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ABSTRACT 


Descriptions  are  given  for  four  programs  which  compute  the  com¬ 
ponents  of  velocity  in  the  wave  train  of  a  unit  source.  The  programs 
give  the  velocity  at  an  arbitrary  position  in  any  direction  and  at 
any  distance  from  the  source.  The  programs  have  been  used  to  com¬ 
pute  a  table  of  components  of  velocity  for  27000  positions. 


ZUSAMMENFASSUNG 


Es  werden  vier  Programme  beschrieben,  die  die  Geschwindigkeits- 
komponenten  im  Wellenzug  einer  Einheitsquelle  berechnen.  Die 
Programme  liefern  die  Geschwindigkeit  fur  einen  beliebigen  Punkt  in 
jeder  Richtung  und  in  jedem  Abstand  von  der  Quelle.  Die  Programme 
sind  benutzt  worden,  urn  eine  Tabelle  fur  die  Geschwindigkeits- 
komponenten  an  27000  Punkten  zu  berechnen. 


RESUME 


On  presente  des  descriptions  pour  quatre  programmes  qui  calculent 
les  composants  de  vitesse  dans  le  train  d'ondes  d'une  source  unitaire. 
Les  programmes  indiquent,  dans  une  position  arbitraire,  la  vitesse 
dans  n'importe  quelle  direction  et  a  n’importe  quelle  distance  de  la 
source.  Les  programmes  ont  ete  utilises  afin  de  constituer  un  tableau 
des  composants  de  vitesse  pour  27.000  positions. 
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FOREWORD 


The  objective  of  this  report  is  the  release  of  enough  information 
so  that  a  programmer  can  utilize  the  existing  programs  for  particular 
applications  without  waiting  for  a  full  documentation  of  the  symbolic 
and  numerical  analysis.  The  work  was  supported  initially  by  the 
Office  of  Naval  Research  under  Project  No.  NR  062-203,  but  is  supported 
currently  as  part  of  the  Foundational  Research  Program  of  the  Naval 
Weapons  Laboratory,  Project  No.  R360FR103/2101/R011010-1.  Assistance 
for  programming  and  checkout  was  contributed  by  W.  H.  Langdon  and 
E.  J.  Hershey.  The  date  of  completion  of  this  report  was  30  June 
1965. 


APPROVED  FOR  RELEASE: 


/s/BERNARD  SMITH 
Technical  Director 
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INTRODUCTION 

The  objective  of  work  on  ship  waves  at  the  Naval  Weapons 
Laboratory  is  the  computation  of  waves  around,  ships  of  finite  breadth 
in  water  of  finite  depth.  Thus  the  scope  of  the  work  is  not  limited 
to  the  far  field  directly  astern  of  the  ship.  A  computation  of  the 
flow  over  the  hull  of  a  ship  is  basic  to  the  computation  of  both  the 
friction  drag  and  the  wave  drag  on  the  ship.  The  computation  would 
help  to  improve  the  efficiency  of  ships.  Unfortunately,  it  has  not 
been  possible  heretofore  to  bring  the  calculator  into  competition 
with  the  towing  tank.  It  costs  more  money  to  set  up  a  model  on  the 
calculator  than  it  costs  to  make  a  model  and  test  it  in  a  towing  tank. 
A  major  effort  has  been  directed  therefore  toward  improvement  of  the 
method  of  computation. 

The  time  consuming  step  in  the  method  is  the  calculation  of 
velocity  in  the  wave  train  which  trails  behind  a  unit  source.  The 
components  of  velocity  are_ derived  from  Havelock's  double  Fourier 
integral  in  wave  number  space.  Following  a  suggestion  from  the  David 
Taylor  Model  Basin,  A.  R.  DiDonato15  has  developed  a  method  of  inte¬ 
gration  in  which  the  path  of  integration  in  the  complex  plane  is  so 
displaced  as  to  avoid  a  singularity  in  the  integrand.  On  the  other 
hand,  radial  integration  through  the  singularity  introduces  the 
complex  exponential  integral  and  simplifies  the  double  integration  to 
single  integration.  Azimuthal  integration  still  must  be  applied" to 
an  integrand  which  is  partly  oscillatory  and  partly  monotonic. 

The  computation  of  special  functions  for  random  arguments  is 
achieved  on  a  digital  calculator  by  reference  to  subroutines.  For 
the  ship  wave  application  the  subroutines  may  require  several  seconds 
for  computation.  The  development  of  subroutines  has  been  continued 
through  many  generations  to  what  is  considered  to  be  the  ultimate 
from  the  standpoint  of  efficiency.  Four  systems  of  subroutines  and 
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programs  are  the  subject  of  the  present  report.  The  programs  are 
coded  in  language  for  the  Naval  Ordnance  Research  Calculator  but  are 
not  available  in  FORTRAN.  The  NORC  does  thirteen  decimal  digit 
arithmetic  at  the  rate  of  15000  operations  per  second  in  response  to 
three  address  instructions.  The  possibility  of  translating  the  pro¬ 
grams  to  other  languages  is  under  consideration. 

Between  the  four  programs  it  is  possible  to  obtain  output  with 
better  than  six  decimal  digit  accuracy,  whereas  the  mathematical 
model  for  ship  waves  is  in  error  insofar  as  it  does  not  meet  exactly 
the  physical  boundary  conditions.  The  components  of  velocity  from 
the  computing  programs  are  to  be  used  in  the  formation  of  a  matrix 
which  may  be  nearly  singular,  and  unless  a  relatively  high  level  of 
accuracy  is  maintained,  the  inversion  of  the  matrix  may  introduce  more 
error  than  the  mathematical  model. 

In  the  first  computing  program,  the  azimuthal  integration  is  com¬ 
pleted  by  the  application  of  an  integration  rule  to  computed  values 
of  the  integrand.  Inasmuch  as  the  angular  variable  of  integration  is 
cyclical,  the  high  accuracy  rule  of  integration  is  the  trapezoidal 
rule  for  equally  spaced  arguments.  This  method  has  been  used  for  sub¬ 
merged  bodies  of  revolution,  but  it  is  impractical  for  surface  ships. 
Thousands  of  intervals  of  integration  are  required,  and  it  may  take  an 
hour  to  evaluate  one  velocity  on  NORC. 

A  breakthrough  has  been  achieved  through  an  integration  by  parts 
which  makes  it  possible  to  integrate  through  many  oscillations  at  a 
time.  It  may  be  recognized  that  the  integrand  contains  two  factors 
one  of  which  is  monotonic  while  the  other  contains  all  of  the  fluctua¬ 
tion  of  the  integrand.  If  the  first  factor  is  approximated  in  terms 
of  the  powers  of  a  parameter  common  to  both  factors,  then  the  approxi¬ 
mation  can  be  integrated  term  by  term  through  the  use  of  recurrence 
equations. 
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In  the  second  computing  program,  the  monotonic  factor  is  approxi¬ 
mated  in  terms  of  positive  integral  powers  of  the  common  parameter. 

This  has  the  advantage  that  only  single  valued  complex  functions  are 
required  for  integration,  and  no  distinction  must  be  made  between 
different  branches  of  the  Riemann  surface.  It  has  the  disadvantage 
that  irrational  functions  must  be  approximated,  and  the  range  of 
approximation  is  limited.  Hundreds  of  intervals  of  approximation  are 
required,  but  it  takes  only  32  seconds  to  evaluate  one  velocity  on  NORC. 

The  classical  approach  to  a  surface  disturbance  by  Kelvin1,2’3  was 
an  application  of  the  principle  of  stationary  phase  to  the  determina¬ 
tion  of  the  first  terms  of  an  asymptotic  expansion.  Considerations  of 
the  possibility  of  development  into  full  expansions  have  been  made  for 
the  line  of  wave  crests  by  Peters12  and  by  Ursell13. 

In  the  present  approach  to  submerged  sources,  the  integrand  is 
expressed  as  the  sum  of  a  monotonic  term  and  an  oscillatory  term.  In 
terms  of  a  common  parameter  the  phase  of  the  oscillatory  term  is  expressed 
as  a  polynomial  of  at  most  the  third  degree,  while  the  amplitude  of  the 
oscillatory  term  is  expressed  as  a  power  series  of  unlimited  degree. 
Symbolic  expressions  for  the  terms  are  very  complicated. 

In  the  third  computing  program,  symbolic  analysis  is  combined  with 
numerical  computation  to  generate  the  series  expansions.  The  asymptotic 
expansion  is  valid  at  large  distances  astern  of  the  source  or  ahead  of 
the  source  but  the  Stokes  phenomenon  makes  trouble  abeam  of  the  source. 

A  representative  computing  time  is  6  seconds  to  evaluate  one  velocity  on 
NORC. 


In  view  of  the  relatively  high  efficiency  of  the  asymptotic  expan¬ 
sion,  a  further  effort  was  made  to  improve  the  integration  by  parts.  The 
convolution  theorem  offers  the  best  promise  of  transformation  of  integrals 
without  loss  of  accuracy,  but  it  leads  to  integrals  of  multiple  valued 
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functions  for  which  new  subroutines  have  had  to  be  developed.  The 
invention  of  algorithms  for  keeping  the  path  of  integration  on  con¬ 
tiguous  branches  of  the  Riemann  surface  required  the  guidance  of 
graphical  analysis,  and  the  invention  of  formulae  for  determining 
ranges  of  approximation  for  balanced  accuracy  required  the  guidance 
of  experimental  computations. 

In  the  fourth  computing  program,  the  monotonic  factor  is  approxi¬ 
mated  in  terms  of  positive  and  negative  half  integral  powers  of  a 
parameter  common  to  both  factors.  When  either  of  two  branch  points  in 
the  complex  plane  is  selected  to  be  the  origin  of  the  parameter,  that 
branch  point  is  replaced  by  a  branch  point  elsewhere  on  the  Riemann 
surface,  and  the  range  of  approximation  can  be  expanded.  Only  four 
intervals  of  approximation  are  required,  but  it  still  takes  9  seconds 
to  evaluate  one  velocity  on  NORC. 

The  computation  of  special  functions  at  greater  speed  can  be 
achieved  through  interpolation  in  a  table  of  sufficient  size.  The  time 
for  random  access  in  the  table  increases  while  the  time  for  interpolation 
decreases  with  increase  in  the  size  of  the  table.  For  the  ship  wave 
application  the  table  could  contain  more  than  a  hundred  thousand  numbers. 
A  table  of  this  size  could  not  all  be  stored  in  the  core  memory  of  many 
calculators,  and  the  access  time  from  magnetic  tape  would  determine  the 
time  of  operation. 

For  studies  in  interpolation  a  table  has  been  prepared  with-.27000 
sets  of  components  of  velocity.  The  -table  is  accessible  to  the  IBM  7090 
series  of  calculators.  The  table  has  been  recorded  in  the  binary  coded 
decimal  mode  on  magnetic  tape  with  one  80  column  card  image  per  record. 
This  table  still  is  not  large  enough  to  provide  a  universal  basis  for 
interpolation,  and  a  further  enlargement  of  the  table  is  under  considera¬ 
tion. 


4 


NWL  REPORT  NO.  198? 


Interpolation  is  especially  difficult  where  both  transverse  and 
divergent  wave  systems  are  present  in  the  wave  train.  Computation 
determines  the  sum  of  the  wave  systems,  but  otherwise  their  ampli¬ 
tudes  and  phases  are  arbitrary.  Separation  of  the  wave  systems  by 
analysis  leads  to  asymptotic  series  which  are  accurate  only  at  large 
distances  from  the  source.  Convergence  of  the  asymptotic  series  can 
be  accelerated  through  the  use  of  converging  factors,  but  experience 
has  indicated  that  converging  factors  themselves  can  be  asymptotic. 
Development  of  rational  functions  in  the  range  of  difficulty  has  been 
unsuccessful  because  zeros  appeared  in  the  denominators  of  the 
rational  functions. 

Documentation  of  the  symbolic  and  numerical  analysis  is  under 
preparation  but  will  require  a  long  time.  In  the  meantime,  it  is 
easy  to  use  any  of  the  programs  for  computations.  Only  the  first  two 
programs  provide  for  water  of  finite  depth,  although  all  four  programs 
provide  for  water  of  infinite  depth.  The  present  report  is  restricted 
to  the  case  of  infinite  depth,  and  a  later  report  will  take  up  the 
case  of  finite  depth. 


ANALYSIS 


Formulae  for  the  velocity  in  the  wave  train  of  a  point  source  have 
been  derived  for  infinite  depth  by  Havelock10  and  for  finite  depth  by 
Lunde11.  A  brief  derivation  is  reproduced  herein  as  part  of  the  docu¬ 
mentation  of  notation. 

Let  the  point  source  be  moving  at  a  depth  h  below  the  surface  with 
a  speed  U.  A  point  in  the  fluid  is  located  at  Cartesian  coordinates 
x,  y,  z  in  a  righthanded  coordinate  system  which  is  moving  with  the 
source.  The  origin  of  coordinates  is  at  the  surface  directly  over  the 
source.  The  x-coordinate  is  the  distance  forward,  the  y-eoordinate  is 
the  distance  to  the  right,  and  the  a- coordinate  is  the  distance  down¬ 
ward. 
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In  accordance  with  the  usual  assumption  of  incompressible,  irro- 
tational  flow,  the  velocity  in  the  fluid  is  expressed  as  the  negative 
gradient  of  a  velocity  potential.  Under  steady  state  conditions  the 
potential  in  the  moving  reference  frame  is  given  by  the  expression 

Ux  +  <p(x.  y,  z )  (l) 


where  Ux  is  the  potential  for  uniform  flow  in  a  direction  opposite 
to  the  motion  of  the  source  and  (p(x,  y,  z)  is  the  potential  for  the 
local  disturbance  in  the  wave  train.  The  computing  routines  give  the 
Cartesian  components  u,  v,  tc  of  velocity  in  accordance  with  the 
equations 


3cp 

u  -  -  — 
dx 


V 


3<p 

dy 


3<p 

10  - - 

dz 


(2) 


Let  the  configuration  of  the  free  surface  be  expressed  by  the 
equation 


2  -  f(x,  y)  =  0 


(3) 


For  steady  state  conditions  the  velocity  at  the  free  surface  is 
tangent  to  the  surface  and  the  potential  obeys  the  boundary  equation 


dx  dx  dy  dy  dz 

The  equation  is  linearized,  through  neglect  of  the  term 

dw  df  dcp  df 

-  -  -|-  -  -  (V  Q 

dx  dx  dy  dy 


to  give  the  boundary  equation 


(4) 


(5) 


U^l  _  ^ 

dx  dz 


(6) 
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Let  the  motion  of  the  fluid,  be  determined  by  the  Bernoulli  equation, 

;{(*+£>'  +  <f?>2+  (f?)2}  +  t  -  s,  =  constant  (7) 
2  l  dx  dy  dz  p  ' 

where  p  is  the  density,  j>  is  the  pressure,  and  g  is  the  acceleration 
of  gravity.  At  the  free  surface  the  pressure  is  constant,  and  the 
equation  is  linearized,  through  neglect  of  the  term 

,  dcp  2  dcp  .2  .  Sop  .2  .  . 

dx  dy  dz 

to  give  the  boundary  equation 

Up  -  gf  =  0  (9) 

dx 

The  free  surface  is  eliminated  from  Equations  (6)  and  (9)  by  dif¬ 
ferentiation  to  give  the  equation 

p2  -  ^  =  0  (10) 

dx2  0  dz 

where  the  critical  wave  number  x0  is  defined  by  the  equation 


0  U2 


Along  the  line  behind  the  source  the  wave  length  X  of  the  transverse 
waves  is  given  by  the  equation 


(12) 


7 


ML  REPORT  NO.  1987 


and  a  Froude  number  based  upon  a  length  l  is  given  by  the  equation 


V  gl 

Analysis  shows  that  the  potential  9  may  be  expressed  as  the  sum  of 
three  potentials  in  accordance  with  the  equation 

<p(x,  y,  z)  =  %.(*,  y,  z)  +  tp2  (x,  y,  z )  +  %(x,  y,  z)  (14) 

where  is  the  potential  of  the  source  in  an  unbounded  fluid,  qp2  is 
the  potential  of  an  image  source  over  the  free  surface,  and  cp3  is  the 
potential  of  the  wave  train. 

The  potential  is  given  by  the  equation 


1 

<Pl  =  - T 

{x2  +  y2  +  ( z-h )2}2 

and  its  derivatives  are  given  by  the  equations 


Scpi 

X 

dx 

{x2 

+  y2 

+ 

(2- 

-*)*)* 

9cpi 

y 

dy 

{x2 

+  y2 

+ 

(*- 

-h)2f 

z 

- 

h 

dz 

{x2 

+  y2 

+ 

(*- 

-h)2)* 

(15) 


(16) 


(17) 


(18) 
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In  the  particular  case  where  the  function  is  given  by  the  equation 


F(x,  y)  - 


2  _L  „2  r 


✓  *  +  y 


(25) 


the  amplitude  is  given  by  the  equation 


-i  271  OS  . 

1  „  „  -tnrcostf) 


i(x,  6)  =  — 2  /  f  e 


ATI 


drd(fi 


o  o 


(26) 


where  r,  cfr  are  polar  coordinates  in  coordinate  space. 


It  may  be  deduced,  with  the  aid  of  the  equations 


14 


1  n 

Jn(z)  =  -  f  cos (z  cos  8)  d6 

n  n 


(27) 


and 


/  J0(z)dz  =  1 


(28) 


that  the  amplitude  is  given  by  the  equation 


i(x,  e)  =  ~ 

'271X 


(29) 
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Generalization  to  a  three  dimensional  solution  of  Laplace's  equation 
leads  to  the  equation 


y,  z) 


1 

■2% 


+7T  “ 

/•  / 


-71  0 


-X 

e 


3 


ft 


+  tx(xcos0  +  vslnS ) 

dydd 


(30) 


for  the  potential  of  the  point  source  in  an  unbounded  fluid.  To 
this  must  be  added  a  term  such  that  the  sum  of  terms  satisfies  the 
boundary  condition.  Differentiation  and  substitution  shows  that  the 
amplitude  of  the  additional  term  is  given  by  the  equation 


a(k,  e)  = 


1  (x0  +  X  COS2  6) 

271 X  (x0  -  X  cos2  d) 


or  after  rearrangement,  by  the  equivalent  equation 


a(x,  e) 


1  +  *o _ 

27tX  7tx(x0  -  X  COS2  6) 


(31) 


(32) 


The  first  term  in  this  expression  for  amplitude  leads  to  the  equation 


.  .1  r°  -x(a  +  ft)  +  ix(xcos0  +  ysln0 )  , 

cp2(x,  y,  z)  = - J  f  e  dy.d6  (33) 

-Tj  o 


for  the  potential  of  the  image  source  over  the  free  surface.  The 
second  term  in  the  expression  for  amplitude  leads  to  an  equation 
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which  can  be  simplified  through  the  substitutions 


(z  +  h)  -  i(x  cos  6  +  y  sin  6) f 


and 


x 


x0 

cos2  6 


u 

(z  +  h)  -  i(x  cos  6  +  y  sin  8) 


(34) 


(35) 


The  function  e  satisfies  both  Laplace's  equation  and  the  boundary- 
equation  for  the  free  surface.  It  is  added  in  just  the  right  amount 
to  make  the  wave  train  trail  behind  the  source  if  the  path  of  inte¬ 
gration  with  respect  to  u  proceeds  along  a  straight  line  in  the 
complex  plane  from  «  =  -  ®  toward  the  origin,  bypasses  the  origin  on  a 
half  circle  of  small  radius  below  the  origin,  and  continues  on  a  straight 
line  to  u  =  5.  The  integrand  is  analytic  and  the  path  of  integration  may 
be  deformed  in  accordance  with  the  Cauchy  theorem  on  analytic  functions. 
The  potential  of  the  wave  train  is  given  by  the  equation 

K  +71  e-8  Sew 

%  (x,  y,  z)  -  —  S  - —a  f  -  dude  (36) 

n  cos2  6  u 


Differentiation  of  this  potential  requires  the  equation 


d. 

db 


—  00 


-8  6  eu 

e  S  -  du 


(37) 
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which  may  be  integrated  by  parts  to  give  the  equation 


d  -6  8  eu 


d  8 


8  „6  eu 


e  f  -  du  =  -  e  "  J  -0  du 


(38) 


The  derivatives  of  <ps  are  given  by  the  equations 


_ 
dx 

_ 
dy 

_ 

9a  Ti 


i  *0 

+71 

r 

e-5 

71 

j 

-71 

cos4  8 

.  x| 

+71 

r 

e-8 

n 

j 

“7C 

cos4  6 

*0  r 

+7C 

‘-6  , 

8  ew 


4—,  /  ~2  du  cos  0  90 


8  eu 


cos4  9 

-7C  -CO 


(39) 


(40) 


(41) 


Only  the  real  parts  of  the  complex  computations  are  reported  by  the 
computing  programs. 


PROGRAMMING 


The  computation  on  NORC  requires  four  magnetic  tapes. 

The  input  is  on  tape  code  01.  The  input  consists  of  four  word 
blocks  which  are  numbered  serially  in  ascending  order  of  block  number. 
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Each  block  of  input  contains  a  set  of  values  for  the  coordinates 


h,  x,  y,  z 


The  blocks  of  input  are  terminated  with  an  end  of  file. 

The  output  is  on  tape  code  02.  The  first  block  of  output  is  num¬ 
bered  9999  and  contains  the  six  quantities 


x0>  d,  n,  0,  0,  0 


where  x0  is  the  critical  wave  number,  d  is  the  depth  of  water,  and  n 
is  the  number  of  intervals  for  trapezoidal  integration.  The  remain¬ 
ing  output  consists  of  seven  word  blocks  which  are  numbered  like  the 
input.  Each  block  of  output  contains  a  set  of  values  of  the 
quantities 


h,  x,  y,  z,  u,  v,  w 


The  blocks  of  output  are  terminated  with  an  end  of  file. 

The  library  is  on  tape  code  09.  The  subroutines  for  the  library 
are  in  Deck  2500. 

The  programs  are  on  tape  code  12.  The  blocks  of  programming  are 
in  Deck  2513,  which  is  arranged  as  shown  on  the  following  page. 
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Block  Number  Memory  Address  Program 


9999 

0002 

- 

0007 

Start  #1 

0001 

3861 

- 

3892 

Assembly 

0002 

0704 

- 

0928 

Computation 

0003 

3952 

— 

4000 

Input  -  Output 

9998 

0002 

- 

0007 

Start  #2 

0001 

3861 

- 

3892 

Assembly 

0002 

0704 

- 

1355 

Computation 

0003 

3952 

- 

4000 

Input  -  Output 

9997 

0002 

- 

0007 

Start  #3 

0001 

3873 

- 

3892 

Assembly 

0002 

0420 

- 

1159 

Computation 

0003 

3937 

- 

4000 

Input  -  Output 

9996 

0002 

- 

0007 

Start  § 4 

0001 

5001 

- 

5221 

Assembly 

0002 

2989 

- 

3991 

Computation 

0003 

4940 

- 

5000 

Input  -  Output. 

End  of  File. 

Each  computing  program  consists  of  a  set  of  three  program  blocks 
which  are  preceded  by  a  starting  block.  To  select  a  program  it  is  only 
necessary  to  key  into  location  0001  an  instruction  which  reads  in  the 
starting  block  and  then  start  the  calculator  with  the  instruction  in 
location  0001.  The  starting  instruction  for  each  of  the  four  programs 
is  as  follows. 

Program  Starting  Instruction 


#1. 

Trapezoidal  Rule 

1294 

0002 

0007 

9999 

#2. 

Integration  by  Parts 

1294 

0002 

0007 

9998 

#3. 

Asymptotic  Series 

1294 

0002 

0007 

9997 

#4. 

Integration  by  Convolution 

1294 

0002 

0007 

9996 
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At  the  start  of  a  run,  the  assembly  routine  reads  into  memory 
addresses  3894  -  4000  an  assembly  subroutine  from  the  library.  The 
assembly  routine  assembles  subroutines  from  the  library,  reads  in  the 
rest  of  the  program,  and  transfers  control  to  the  input-output 
routine.  The  space  required  for  subroutines  by  the  four  programs  is 
as  follows. 


Program 


Memory  Addresses 


§1.  Trapezoidal  Rule 
§2.  Integration  by  Parts 
#3.  Asymptotic  Series 
#4.  Integration  by  Convolution 


0008  -  0703 
0008  -  0703 
0008  -  0418 
0008  -  2988 


During  the  progress  of  computation,  numbers  are  accumulated  and  stored 
temporarily  in  common  storage  in  the  memory.  The  space  required  for 
common  storage  by  the  four  programs  is  as  follows. 


Program 

Memory  Addresses 

#1- 

Trapezoidal  Rule 

0929  - 

3951 

#2. 

Integration  by  Parts 

1356  - 

3951 

#3. 

Asymptotic  Series 

1248  - 

1999 

#4. 

Integration  by  Convolution 

3998  - 

4936 

The  calculator  starts,  restarts,  or  stops  under  the  control  of 
switches.  Computation  starts  at  the  beginning  of  the  output  tape  when 
Switch  74  is  turned  off,  but  restarts  after  the  last  block  written  on 
the  output  tape  when  Switch  74  is  on  transfer.  The  computation  continues 
as  long  as  Switch  77  is  on  transfer,  but  makes  a  standard  stop  after  it 
has  completed  a  block  of  output  whenever  Switch  77  is  turned  off.  The 
programs  all  operate  only  in  10K  storage. 
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TABULATION 

The  results  of  computation  on  NORC  are  recorded  card  by  card  on 
magnetic  tape.  On  each  card  the  first  11  columns  are  blank,  columns 
12  to  75  contain  four  16  digit  words,  and  the  last  5  columns  are 
blank. 

Each  number  word  contains  16  digits  such  that  the  first  two 
digits  give  the  exponent  to  base  10  (with  negative  exponents  indi¬ 
cated  by  9's  complements),  the  third  digit  gives  the  sign  (with  0 
for  +  and  1  for  -),  while  the  remaining  thirteen  digits  give  the 
decimal  coefficient,  with  the  decimal  point  after  the  first  of  the 
thirteen  digits.  The  data  are  organized  into  sets,  with  values  of 
h,  x,  y,  z,  u,  v,  w  in  each  set. 

The  sets  are  packed  into  blocks  such  that  the  first  and  last 
columns  contain  an  extra  punch  in  row  11  and  the  first  and  last  words 
are  beginning  of  block  or  end  of  block  words.  Each  BOB  or  EOB  word 
contains  four  fields  of  4  digits  each  such  that  the  PQ  field  is  0290 
or  0291,  the  R  field  is  the  first  memory  address  of  information  in 
the  block;,  the  S  field  is  the  last  memory  address  of  information  in 
the  block,  and  the  T  field  is  the  block  number.  Each  block  contains 
90  sets  of  data  or  630  number  words  for  90  distances.  There  are  two 
groups  of  distances,  for  r  <  1  and  for  r  >  1,  within  each  of  which 
the  spacing  is  similar  to  the  spacing  for  Chebyshev  approximation. 

Each  block  gives  a  different  azimuth  out  of  a  set  of  31  azimuths. 

There  are  three  groups  of  azimuth,  for  n  ^  4>  $  sin-1  j,  for 
sin-1  >  4>  >  |n,  and  for  \ %  >  4>  >  0,  within  each  of  which  the  spacing 
is  the  same  as  the  spacing  for  Lobatto  approximation.  The  blocks  of 
data  are  organized  into  files. 

Each  file  starts  with  a  block  numbered  9999  which  states  that 
x0  =.  1  and  d  =  oo.  The  blocks  of  data  are  numbered  0001  to  0031. 
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The  file  ends  with  an  EOF  word  which  gives  the  file  number.  Each  file 
gives  a  different  depth  out  of  a  set  of  10  depths,  for  g  h  <  16, 
where  each  consecutive  depth  is  twice  the  previous  depth. 

The  inequality  of  spacing  was  intended  to  provide  a  more  accurate 
basis  for  interpolation  than  is  possible  with  equality  of  spacing. 

DISCUSSION 

In  the  classical  analysis  of  wave  resistance  by  Michell4  the  poten¬ 
tial  was  expressed  as  the  sum  of  three  terms  of  which  two  cancelled  out 
in  the  integration  over  the  surface  while  the  third  term  alone  con¬ 
tributed  to  the  wave  resistance.  Analogous  cancellations  arise  when 
the  present  functions  are  applied  to  the  limiting  case  of  a  thin  ship. 

Let  da  be  defined  to  be  the  projection  on  a  cross  sectional  plane 
of  one  surface  element  on  the  side  of  the  ship.  In  the  limiting  case 
of  a  thin  ship  the  Bernoulli  pressure  on  this  element  would  contribute 

-  p Up  da  (42) 


to  the  wave  resistance.  •  The  total  potential  is  the  sum  of  contribu¬ 
tions  from  other  surface  elements.  Let  da'  be  defined  to  be  the  pro¬ 
jection  on  a  cross  sectional  plane  of  another  surface  element.  The 
flux  through  the  projection  da’  would  be  Uda ’  for  uniform  flow.  This 
flux  is  cancelled  at  the  surface  by  a  local  source  of  strength 


0  ,  \ 
—  da'  (43) 

471 
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on  the  median  plane.  The  total  wave  resistance  R  is  given  by  the 
double  integral 


R 


V) 


p  u2 

An 


dado' 


(44) 


where  the  integration  with  respect  to  dado'  is  taken  over  both  port 
and  starboard  halves  of  the  hull. 

For  every  pair  of  elements  which  constitute  a  pressure  point  and 
a  source  point  there  is  another  pair  of  elements  with  the  pressure 
point  and  the  source  point  interchanged,  and  for  which  the  coordinate 
x  in  the  expression  for  cp  is  reversed  in  sign.  Any  term  in  <p  for 
which  dy/dx  is  odd  with  respect  to  x  will  contribute  nothing  to  the 
wave  resistance.  This  disposes  of  the  terms  and  <jp2.  In  the  case 
of  cp3,  a  reversal  of  x  with  y  =  0  replaces  the  path  of  integration 
with  respect  to  a  by  a  complex  conjugate  path  everywhere  except  on 
the  half  circle  which  bypasses  the  origin.  The  only  contribution  to 
the  wave  resistance  comes  from  the  integration  at  the  origin.  The 
surviving  term  is  retained  if  the  potentials  are  replaced  in  accord¬ 
ance  with  the  substitution 


dx 


Hq  / 


+7t  e-5 


-71 


cos®  8 


d6 


(45) 


A  change  in  variable  of  integration  in  accordance  with  the  equation 


cos  8  - 


1 

X 


(46) 
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leads  to  the  equation 


d6 


”  -\2x0(z+h)  ,  X2dX 

-  J  e  cos  \x0x  — - 

i  /  X2  -  1 


(47) 


where  the  integral  on  the  right  may  be  recognized  as  Michell's  function 
as  defined  by  Birkhoff  and  Kotik6.  ?»7ith  due  recognition  that  the 
Michell  integral  is  expressed  by  an  integration  once  only  over  the 
median  plane,  it  may  be  verified  that  the  present  formulae  are  con¬ 
sistent  with  Michell's  integral  in  the  limiting  case  of  a  thin  ship. 

It  might  be  mentioned  that  a  computing  program  for  generating 
Michell's  function  along  with  other  oscillatory  parts  of  the  velocity 
has  been  constructed  incidental  to  efforts  to  achieve  a  successful 
method  for  interpolation. 

The  Michell  function  is  a  special  case  of  a  set  of  functions  which 
have  been  investigated  by  Wilson7.  The  functions  (cc,  (3)  and 
(a,  (3)  are  the  real  and  imaginary  parts  of  a  complex  function 
E^  (a,  |3)  as  stated  by  the  equation 

E(#)  (a,  j3)  =  4n)  (<*-  P)  +  (a-  P)  (48) 

while  the  complex  function  is  defined  by  the  equation 


n 


E(*}  (a,  p)  =  / 


-a2  sec 2  8 


iBsec#  _ 
sec 


e  de 


(49) 


The  functions  of  different  orders  have  convergent  and  asymptotic  expan¬ 
sions  and  they  are  connected  by  recurrence  equations.  Values  of 
E^  (a,  {3)  for  n  -  -  1,  0,  1,  2,  3  have  been  tabulated  in  reference  (8). 
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The  Michell  function  is  the  real  part  of  the  third  order  Wilson 
function.  A  quantitative  comparison  with  functions  in  the  present 
report  is  possible  only  for  (3  =  0.  The  substitution 


sec  6  -  -  (l  +  t) 
2 


(50) 


and  evaluation  in  terms  of  Macdonald  functions  gives  the  equation 


8 ,  14 


(a,  0)  =  -  e  2  f  e  2 


1  -#2,”  ~%2t  dt  1 


V  t2  -1 


=  -  e  2  K0(- 1  a2)  (51) 


whence  differentiation  with  respect  to  a2 gives  the  equation 


B,  14 


4S)  (a,  0)  =  i  «■“  (J-0(i«2>  +  Mi<x2» 


(52) 


At  x  -  y  =  0  the  path  of  integration  with  respect  to  u  lies  along  the 
real  axis  except  at  the  half  circle  where  it  bypasses  the  origin. 
Integration  with  respect  to  u  and  the  substitutions 

a2  =  >c0(z  +  h)  P  =  k0x  (53) 

lead  to  the  equation 


d% 

dx 


4x2  (a,  0) 


(x  =  y  =  0)  (54) 


This  equation  is  fulfilled  to  ten  digit  accuracy  by  the  numerical 
integrations. 
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At  the  ONR  seminar  on  wave  resistance,  Inui8  has  described  tables 
from  which  the  full  value  of  dcpB/dx  may  be  computed  with  both  monotonic 
and  oscillatory  parts  included,  but  only  for  the  direction  astern  of 
the  source. 


CONCLUSION 

Any  system  which  is  used  to  compute  velocity  for  a  three  dimen¬ 
sional  position  must  respond  to  a  random  reference.  Subroutines  can 
be  utilized  for  an  efficient  evaluation,  but  interpolation  may  require 
too  bulky  a  table-  The  fastest  subroutines  require  6  to  9  seconds  per 
evaluation  on  NORC. 
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